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i n t roduct ion 

The barotropic finite element code described by Staniforth 
and Mitchell in Ref. 1 has been the point of departure for exten- 
sive efforts at the Naval Postgraduate School to develop improved 
techniques for numerical weather prediction. Two earlier reports 

(Refs. 2 & 3) described applications of tensor product techniques 
\ 

to extensions of the code developed by Hinsman (Ref. 4). The 
present report is concerned with proposed alterations of the 
Stani f or th-Mi tche 1 I code for the purpose of improving computa- 
tional efficiency. 

Three separate contributions are reported here. First of 
these is a scheme for the direct solution of the Helmholtz equa- 
tion as a substitute for the Concus and Golub iterative algorithm 
(Ref. 5) used in the Stani for th-Mi tche 1 1 code. The second item 
is an improved solution of the generalized eigenvalue problem 
that arises in transforming the two-dimensional Poisson and Helm- 
holtz problems into a series of one-dimensional problems. This 
solution takes advantage of the special form of the matrices 
involved (symmetric, tridiagonal) and avoids matrix inversion. 
The third item deals with the special form of the generalized 
eigenvalue problem that arises when there is a periodic boundary 
condition in the east-west direction. In this instance the 
symmetric, formerly tridiagonal, matrices remain symmetric, but 
have nonzero elements in the upper right-hand and lower left-hand 
corners. A separate solution algorithm is required for determin- 
ation of the eigenvalues and also for the eigenvectors. 
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D i rect So 1 u t i on of the He 1 mho 1 t z Equal i on 



The Helmholtz equation may be written in the form 

delsqw-hw=f <1> 

where delsq is the two-dimensional Laplace operator, h is a func- 
tion of y, f is a function of x and y, and w is the dependent 
variable. For the Finite Element (FE) discretization, consider 
the grid of Fig. 1. There are e east-west grid lines and n north- 




e rows 



n colurnns 



(Mote the periodic 
boundary condition 
in the E-W direction. ) 



Fig. 1. Node numbering and spacing. 



south grid lines, defining ne intersections (nodes). Node num- 
bering is from west to east, beginning at the southwest corner. 
The positive x direction is eastward and the positive y direction 
is northward. In the Galerkin FE discretization process the pro- 
duct h w is treated as a single entity and the tensor product 
concept is employed in writing the matrix form of the equation. 
The r esu It is 

MX W SY + SX W MY - MX W H MY = V <2> 

where MX and SX are n x n symmetric "mass" and "stiffness" mat- 
rices for the X direction, MY and SY are the corresponding e x e 
matrices for the y direction, W i s an n x e matrix of nodal 
values of the dependent variable, H is diagonal e x e, and V is 
n X e. MX and SX depend only on the node spacing ai and MY and 
SY depend only on the node spacing bj . Defining equations for 
these matrices are given in Appendix A. Note that matrices SX 
and SY are the negatives of standard FE stiffness matrices. The 
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columns of W contain nodal values of w from the west-east rows, 



beginning with the most southerly. If the nodal values of f are 
similarly arranged in an n x e matrix F, then V is the product 
MX F MY. The nonzero entries in the diagonal matrix H are the 

values of h for the successive rows of nodes, beginning with the 
most southerly. 

It is advantageous for the succeeding manipulations to 
transpose the terms of <2> to obtain 

SY WT MX + MY .WT SX - MY H WT MX = VT <3> 

where WT and VT are the respective transposes of W and V. Note 
that all of the other matrices are symmetric. 

Consider the generalized eigenvalue problem 

SX pt = 1, MX p, <4> 

where p, is the ith eigenvector (n x 1) and I, is the correspond- 
ing eigenvalue. Assembling the eigenvectors in an n x n modal 
matrix P and the eigenvalues in an n x n (diagonal) spectral 
matrix L, the result may be exhibited as 

SX P = MX P L <5> 

It is advantageous to require that the eigenvectors be normalized 
so that PT MX P = I, where PT is the transpose of P and I is the 
nth order identity matrix. If <3> is postmu 1 t i p 1 i ed by P and <5> 
is used to replace SX P in the second terra, the result is 

SY WT MX P + MY WT MX P L - MY H WT MX P = VT P <6> 

Let Q = WT MX P and substitute in <6> to get 

SY’ Q + MY Q L = U <7> 

where SY’ = SY - MY H and U = VT P. Observe that SY’ is not 
symmetric, but is tridiagonal. Now the ith column q, of Q may be 
found by solving 
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(SY’ + \, MY) qj = <8> 

where Ui is the ith column of U. The solutions for the n qi may 
be assembled into Q. It is easy to show that 

W = P QT <9> 

where QT is the transpose of Q. This completes the solution. 

Genera 1 i zed Ei genva 1 ue Problem for Symmet r i c Tr idi agona 1 Matr ices 

Solution of the generalized eigenvalue problem <‘4> for the 
symmetric tridiagonal • matrices SX and MX is based on the Sturm 
sequence property and the method of bisection. The theory is 
given by Wilkinson in Ref. 6. Using this theory and generalizing 
an ALGOL program constructed for the standard eigenvalue problem 
(Ref. 7), SUBROUTINE TR 1 E I G has been written to find the eigen- 
va 1 ues . 

For ’any eigenvalue, the corresponding eigenvector can be 
found by rewriting <A> in the form 

(SX - li MX) qi = 0 <10> 

If the first component of qi is chosen as unity, the first scalar 
equation of <10> determines the second component and each suc- 
ceeding scalar equation determines an additional component. The 
vector thus found may then be normalized. SUBROUTINE EIGVEC per- 
forms these operations. 

Once eigenvectors have been determined, SUBROUTINE RAYLEE 
may be used to refine the eigenvalues by evaluating the Rayleigh 
quotient for each eigenvector. Following this an additional call 
to SUBROUTINE EIGVEC enhances the accuracy of the eigenvectors. 
Experience to date indicates that further iterations are super- 
f 1 uous . 
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Genera 1 i zed Ei genva 1 ue Prob 1 em with Per iodic Boundary Cond i t i ons 



When there is a periodic boundary condition in the x direc- 
tion, the tridiagonal form of matrices SX and MX is modified by 
the presence of nonzero elements in the upper right-hand and 
lower left-hand corners. The matrices remain symmetric, but the 
Sturm sequence and bisection procedure described above requires 
modification. Although the underlying strategy is unchanged, 
different software is required. In this instance, the eigen- 
values are found by caMing SUBROUTINE PEREIG. To determine the 
eigenvectors SUBROUTINE EIGVCP is called. A new problem arises 
here, because, if the grid spacing is uniform in the x direction, 
there will be double eigenvalues. Specifically, if n, the number 
of subdivisions in the x direction, is odd, there will be a 
single zero eigenvalue and all of the others will be double. If 
n is even, the eigenvalue of largest abso 1 ute' va 1 ue will also be 
single. Corresponding to double eigenvalues, the eigenvectors 
are not unique and special procedures must be used to construct 
vector pairs having the required orthogonality. SUBROUTINE EIG- 
VCP tests for repeated eigenvalues and constructs the appropriate 
orthogonal pairs of vectors as required. For this periodic 
boundary condition, eigenvalue refinement via the Rayleigh quo- 
tient is effected by calling SUBROUTINE RAYLYP. An additional 
call to SUBROUTINE EIGVCP provides better eigenvectors. Further 
iteration is believed superfluous. 
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Remarks 



Two additional subroutines are needed to utilize the pro- 
posed alterations to the Stani f or th-Mi tche 1 1 program. First of 
these is SUBROUTINE SETUP, a replacement for SUBROUTINE EBVSET. 
The second is SUBROUTINE SLVBVP, a replacement for SUBROUTINE 
EBVP2D. All of the new subroutines are listed in Appendix A. 

The routines described herein have been tested and found to 
perform satisfactorily on a 12 x 12 grid. The parameter RF in 
TRIEIG and PEREIG (currently set at 1.0 E-07) may need to be 
adjusted for larger grids. In SUBROUTINE EIGVCP there is an 
additional parameter, SEP, which also may need adjustment for 
larger grids. It is currently set at 4.0 E-04. 

There is a fundamental problem arising from the use of single 
precision arithmetic in solutions to the Helmholtz equation. 
This results from the fact that the njean geopotential height is 
approximately three orders of magnitude greater than the ampli- 
tude of a representative disturbance. Comparisons between single 
and double precision solutions show clearly that round-off ad- 
versely affects the former. 

Operation counts have been conducted to evaluate the potent- 
ial saving resulting from substituting the Helmholtz direct sol- 
ver for the Concus-Golub iterative scheme. For a 13 x 13 grid 
and assuming only one cycle of iteration, there is a saving of 10 
percent for each time step. For a 90 x 90 grid, the saving is 
reduced to 5 percent. If additional iterations are required for 
the Concus-Golub solution, these savings would, of course, be 
greater. The savings resulting from the proposed eigenvalue sol- 
utions have not been evaluated, but are expected to be substan- 
tial . 
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APPENDIX A. 



MASS” and "STIFFNESS” MATRICES 



MX = i 
(n»4) 



MY = 1 
6 

( e = 3) 



2(ai»+a 1 ) 

a 1 
0 

a 4 

2b‘i 

bi 

0 



a 1 

2(ai+a2) 

az 

0 

bi 

■ 2 ( b 1 +b 2 ) 
b 2 



0 

az 

2 { a 2 +a 3 ) 
a 3 

0 

b2 
2 b 2 



a„ 

0 

a 3 

2 ( a 3 +a 1 * )_ 



SX = 
( n = 4 ) 



i -i 1 

a 4. 3 1 a 1 



a i, 



i .1 .1 i 

ai aia 2 ^ z 



0 



i .i .i i 

d z 32^3 ^3 



i 

a 4 



0 



1 .i _i 

a 3 a 3 a 4 



SY = 
(e = 3) 



i .i .1 i 

bi bib2 bz 



0 



bz 




Notes : 

!• Matrices SX and SY are negatives of those called stiff- 
ness matrices in standard FE usage. 

2. Forms given for MX and SX are for a periodic east-west 
boundary condition. For a Neumann boundary condition, 
the terms containing a, would be omitted. 
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APPENDIX B. FORTRAN LISTINGS 



SUBROUTINE SETUP ( E I GVEC , E I GVAL , B I GE , B I GC , B I GA , S , HX, HY, HELM, DIR, 

1 FOURTH, NI , NJ , CON) 

» SUBROUTINE SETUP IS DESIGNED TO REPLACE EBVSET WHEN THE PERIODIC » 

If BOUNDARY CONDITION IS IMPOSED ON EASTERN AND WESTERN BOUNDARIES » 

C 



c 

r* 


AUTHOR 


: R. E. 


NEWTON, SUMMER 1986 


c 

c 


ARGUMENTS 
OUT - 


EIGVEC 


- MATRIX OF EIGENVECTORS, NI X NI, IN X DIRECTION 


c 




EIGVAL 


- VECTOR OF EIGENVALUES, NI, IN NONDEISCEND I NG ORDER 


c 




B IGE 


- DIAGONAL FACTOR OF COEFFICIENT MATRIX FOR Y 


c 

c 




B IGC 


DIRECTION 

- SUPERDIAGONAL OF UPPER UNIT TRIANGULAR FACTOR OF 


c 

c 




BIGA 


COEFFICIENT MATRIX FOR Y DIRECTION 
- SUBDIAGONAL OF LOWER UNIT TRIANGULAR FACTOR OF 


c 

c 


IN - 


S 


COEFFICIENT MATRIX FOR Y DIRECTION 
- SQUARE OF MAP FACTOR, N I X N J 


c 




HX 


- NODE SPACING IN X DIRECTION 


c 




HY 


- NODE SPACING IN Y DIRECTION 


c 




HELM 


- LOGICAL SWITCH - TRUE FOR HELMHOLTZ PROBLEM 


c 




DIR 


- LOGICAL SWITCH - TRUE FOR DIRICHLET B.C. ON 


c 

c 




FOURTH 


POISSON PROBLEM 

- LOGICAL SWITCH - TRUE FOR FOURTH ORDER SOLUTION OF 


c 


• 




POISSON PROBLEM 


c 




NI 


- X-DIMENSION 


c 

r* 




NJ 


- Y-DIMENSION 


U 

c 

c- 


NOTE: 


N IS 


MAX(NI ,NJ) 



PARAMETER(N=13) 

REAL AS(N) , BS(N) , CS(N) , S(NI ,NJ ) , HX(NI ) ,HY(NJ) , 

1 AR(N) , BR(N) , CR(N) , BIGA(NI , NJ ) , BIGCCNI , NJ ) , B IGE(NI , NJ ) , WU(N) , 

2 EIGVEC(NI , NI ) , EIGVAL(NI ) , AM ( N ) , BM ( N ) , CM ( N ) , CON 
LOGICAL HELM, DIR, FOURTH 

CF = 2. 

IF(FOURTH) CF = 5. 

C SET UP "MASS" MATRIX FOR X DIRECTION (SUBDIAGONAL AM, DIAGONAL BM, 

C SUPERDIAGONAL CM). PERIODIC BOUNDARY CONDITION ASSUMED. IF NOT 

C PERIODIC, CALL SETABC INSTEAD. 

CALL SETABX(AM, BM, CM, HX, CF, NI ) 

C SET UP "STIFFNESS" MATRIX FOR X DIRECTION (SUBDIAGONAL AS, DIAGONAL 
C BS, SUPERDIAGONAL CS). PERIODIC BOUNDARY CONDITION ASSUMED. IF 
C NOT PERIODIC, CALL SETD2 INSTEAD. 

CALL SETD2N(AS,BS,CS, 1. ,HX,NI ) 

C SOLVE THE GENERALIZED EIGENVALUE PROBLEM: K X = Z M X, WHERE K AND 
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C M ARE STIFFNESS AND MASS MATRICES, RESPECTIVELY, AND X IS THE EIGEN- 
C VECTOR CORRESPONDING TO EIGENVALUE Z. PERIODIC BOUNDARY CONDITIONS 
C ASSUMED. IF NOT PERIODIC, CALL TRIEIG INSTEAD. 

CALL PEREIG(EIGVEC, EIGVAL, AS, BS, BM, AM, WU, N I ) 

C SET UP "MASS" MATRIX (SUBDIAGONAL AM, DIAGONAL BM, SUPERDIAGONAL CM: 
C AND "STIFFNESS" MATRIX (SUBDIAGONAL AS, DIAGONAL BS, SUPERDIAGONAL 
C CS) FOR THE Y DIRECTION. 

CALL SETABC( AM, BM, CM, HY, CF, NJ ) 

CALL SETD2(AS,BS,CS, 1. ,HY,NJ) 

I F( . NOT. HELM)GO TO 11 



C FOR THE HELMHOLTZ PROBLEM THE STIFFNESS MATRIX IS MODIFIED BY SUB- 
C TRACT I NG THE CORRESPONDING ELEMENT OF THE MASS MATRIX DIVIDED BY 
C (GZMN TIMES DELT»»2 TIMES THE RELEVANT SQUARE OF THE MAP FACTOR). 



AS( 1 ) 
BS(1) 
CS( 1) 
NJM = 
DO 10 



= 0 . 

= BS(1 ) 
= CS(1) 
NJ - 1 
J=2, NJM 



BM(1)«C0N/S(1, 1) 
CM(1)»C0N/S(1,2) 





AS(J) = AS(J) 


- AM( J)»C0N/S(1, J-1 




BS(J) = BS(J) 


- BM(J)»C0N/S(1, J) 




CS(J) = CS(J) 


- CM(J)«C0N/S(1,J+1 
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CONTINUE 








AS(NJ) = AS(NJ) 


- AM(NJ)*C0N/S(1,NJM 




BS(NJ) = BS(NJ) 


- BM(NJ)«C0N/S(1,NJ) 




CS(NJ) = 0. 






c 


LOOP OVER EIGENVALUES TO 


CONSTRUCT BIGA 
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DO 20 I = 1, NI 








EIG = EIGVALC I ) 






DO 12 J = 1, 


NJ 






AR(J) = 


AS(J) 


+ EIG»AM(J) 




BR(J) = 


BS(J) 


+ EIG»BM(J) 




CR(J) = 


CS( J ) 


+ EIG«CM(J) 


12 


CONTINUE 







BIGC, AND BIGE 



I F( . NOT. HELM. AND. . NOT. DIR. AND. I . EQ. NI )CR( 1 ) = 0. 



C FACTOR COEFFICIENT MATRIX 

CALL SETTRI (BR, CR, AR, AR, BR, CR, NJ ) 

C STORE RESULTS IN BIGE, BIGC AND BIGA. NOTE SIGN CHANGES FOR BIGC 
C AND BIGA. 



16 

20 



DO 16 J = 1, NJ 
BIGE(I,J) = 
B IGC( I , J ) = 
BIGA( I , J) = 
CONTINUE 
CONTINUE 
RETURN 
END 



BR( J) 
-CR( J ) 
-AR( J ) 
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cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE PEREIG( V, X,B,C,D,E,WU,N) 



* 

* 

* 

» 

» 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 



c 



SUBROUTINE PEREIG USES THE METHOD OF BISECTION TO FIND EIGENVALUES * 
FOR THE GENERALIZED EIGENVALUE PROBLEM INVOLVING SYMMETRIC (ALMOST) * 
TRIDIAGONAL MATRICES WITH NONZERO ELEMENTS IN ’CORNERS’ (PERIODIC » 
BOUNDARY CONDITION). IT IS ADAPTED FROM AN ALGOL PROGRAM NAMED » 
’BISECT’ WRITTEN BY BARTH, MARTIN, AND WILKINSON (NUM. MATH. 9, 386-* 
393 (1967)). ’BISECT’ FINDS EIGENVALUES (STANDARD EIGENVALUE PROB- * 
LEM) FOR A SYMMETRIC TRIDIAGONAL MATRIX. THE MATRICES FOR PEREIG * 
ARE INPUT AS VECTORS - C(N) AS THE DIAGONAL AND B(N) AS THE SUB- » 
DIAGONAL OF THE ’’STIFFNESS” MATRIX AND D(N) AS THE DIAGONAL AND E(N)* 
AS THE SUBDIAGONAL OF THE "MASS" MATRIX. B(l) AND Ed) ARE THE * 
"CORNER” ELEMENTS OF THE RESPECTIVE MATRICES. SUBROUTINE EIGVCP » 
FINDS THE EIGENVECTORS AND NORMALIZES THEM WITH RESPECT TO THE MASS * 
MATRIX. SUBROUTINE RAYLYP USES THR RAYLEIGH QUOTIENT TO OBTAIN IM- * 
PROVED EIGENVALUES USING THESE EIGENVECTORS. A SECOND CALL TO * 
TO EIGVCP EFFECTS A CORRESPONDING IMPROVEMENT IN THE EIGENVECTORS. * 



AUTHOR: R. E. NEWTON, SUMMER 1986 



ARGUMENTS 

OUT - V 

X 

IN - C 
. B 
D 
E 

WU 

N 



MATRIX OF EIGENVECTORS, N X N, 

WITH RESPECT TO "MASS" MATRIX 

VECTOR OF EIGENVALUES IN NONDESCENDING 

DIAGONAL OF "STIFFNESS" MATRIX 

OF "STIFFNESS" MATRIX 
"MASS" MATRIX 
OF "MASS" MATRIX 



NORMALIZED WITH 



ORDER 



SUBDIAGONAL 
DIAGONAL OF 
SUBDIAGONAL 
WORK VECTOR 
VECTOR SIZE 



(= NI ) 



NOTE 



MATRICES ARE FOR X-DIRECTION WITH PERIODIC BOUNDARY 
CONDITION. NODE SPACING IS HX. 



INTEGER N, Z, I , A, K, Nl 

REAL C(N),B(N),X(N),WU(N), EPSl , EPS2,RF, FI , XM I N , XMAX, 

1 Xl,XU,XO,D(N),E(N) ,DC,DB,DD,DE,TMAX,TMIN, V(N,N) , 

2 Q, R, S, Ql, R1 , R2, DR, QTEMP, ROLD, SW 
DATA EPSl, RF/0. , 1. E-7/ 

N1=0 
Z = 0 

CALCULATION OF XMAX, XM I N 



DC=C (N) 

DB=ABS(B(N) ) 

DD=D(N) 

DE=E(N) 

XMAX=(DC+DB) / (DD+DE) 

XMIN=(DC-DB)/ (DD-DE) 

C EIGENVALUES ASSUMED NEGATIVE. IF EIGENVALUES ARE ALL POSITIVE, 
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C REPLACE THE TWO PRECEDING LINES WITH THE FOLLOWING TWO LINES. 

» XMAX= (DC+DB) / (DD-DE) 

XMIN= (DC-DB)/(DD + DE) 

DO 2 I=N-1, 1, -1 
DC=C( I ) 

DB = ABS(B( I ) )+ABS(B( I +1) ) 

DD=D( I ) 

DE=E( I ) +E( I +1 ) 

TMAX= (DC+DB) /(DD+DE) 

TMIN= (DC-DB) / (DD-DE) 

C EIGENVALUES ASSUMED NEGATIVE. IF EIGENVALUES ARE ALL POSITIVE, 
C REPLACE THE PRECEDING TWO LINES WITH THE FOLLOWING TWO LINES. 

» TMAX=(DC+DB)/(DD-DE) 

* TMIN=(DC-DB)/(DD+DE) 

I F(TMAX. GT. XMAX)XMAX=TMAX 
IF(TMIN. LT.XMIN)XMIN=TMIN 
2 CONTINUE 

C SET EPS2 

IF(XMIN+XMAX. GT. 0. )THEN 
EPS2=RF»XMAX 

ELSE 

EPS2=RF» (-XMIN) 

END IF 

IF(EPS1. LE. 0. )EPS1=EPS2 
EPS2=0.5»EPSl+7. »EPS2 

C INNER BLOCK 

XO=XMAX 
DO 4 1=1, N 

■ X( I ) =XMAX 
WU( I )=XMIN 
4 CONTINUE 

C LOOP FOR K-TH EIGENVALUE 

DO 100 K=N, 1,-1 
XU=XMIN 
DO 6 I=K, 1, -1 

IF(XU.LT.WU( I ) )THEN 
XU=WU( I ) 

GO TO 10 
END IF 

6 CONTINUE 

10 IF(XO. GT. X(K) )XO=X(K) 

20 Xl= (XU+XO) /2. DO 

Z = Z + 1 

C STURM’S SEQUENCE 



A = 0 

Q=C( 1 ) -X1»D( 1 ) 
Q1=C(2) -X1*D(2) 
R=B(1)-X1»E(1) 
S=C(N) -X1*D(N) 
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R1=B(2)-X1»E(2) 

R2=0. 

DR = 0. 

DO 30 l=2,N-2 

IF(Q. EQ. 0. )THEN 
SW=1. 

IF(Ql+2. *R1 . EQ. 0. )SW=-1 . 

Q=Q1+SW#2. »R1 
R1=R1+SW*Q1 

R2=SW»(B( I +1 ) -X1*E( I +1 ) ) 

END IF 

IF(Q. LT. 0. )A=A+1 

S=S-R»R/Q 

ROLD=R 

R=DR-R1»R/Q 

QTEMP=Q1-R1»R1/Q 

R1=B( I+1)-X1*E( I+1)-R1»R2/Q 

Q1=C( I +1 )-Xl»D( I +1 ) -R2»R2/Q 

DR=-R0LD»R2/Q 

R2=0. DO 

Q=QTEMP 

30 CONTINUE 

IF(Q. EQ. 0. DO) THEN 
SW=1 .DO 

IF(Q1+2.D0»R1.EQ.0.D0)SW=-1.D0 

Q=Q1+2.D0»SW«R1 

R1=R1+SW»Q1 

R=R+SW* (B(N) -X1»E(N) ) 

IF(Q.LT.0.D0)A=A+1 

ELSE 

IF(Q. LT. 0. ) A=A+1 
END IF 
S = S-R^(R/Q 

R=B(N) -X1#E(N)-R«R1/Q 
Q1=Q1-R1»R1/Q 

IF(Ql.EQ.O. AND.R.NE.O. ) THEN 
A = A+1 

ELSE 

IF(Q1.LT. 0. )A=A+1 
IF(S-R*R/Ql.LT.O. )A=A+1 
END IF 

IF(A.LT.K)THEN 

IFCA.LT. DTHEN 
WU(1)=X1 
XU = X1 

ELSE 

WU(A+1 ) =X1 
XU=X1 

IF(X(A) . GT. XI )X(A)=X1 
END IF 

ELSE 

X0=X1 
END IF 

IF(X0-XU.GT.2. »RF»(ABS(XU)+ABS(XO) )+EPSl)GQ TO 20 
X(K)=(X0+XU)/2. 

100 CONTINUE 

CALL EIGVCP(V, B,C, D, E,X, N) 



13 



CALL RAYLYPCX, V,B,C, WU,N) 

CALL E1GVCP(V,B,C,D,E,X,N) 

RETURN 

END 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx> 

SUBROUTINE E I GVCP ( V, B , C , D , E, X , N ) 

» SUBROUTINE EIGVCP FINDS EIGENVECTORS BY DIRECT SOLUTION OF THE GOV- 

* ERNING LINEAR ALGEBRAIC EQUATIONS. FOR ANY PAIR OF EQUAL EIGEN- 

» VALUES, THE CORRESPONDING VECTORS ARE MADE ORTHOGONAL WITH RESPECT 

* TO THE "MASS” MATRIX. ALL VECTORS ARE NORMALIZED WITH RESPECT TO 

* THE "MASS” MATRIX. 

C 

C AUTHOR - 
C 

C ARGUMENTS 
C OUT - 

C I N - 

C 
C 
C 
C 
C 
C 

c 



R. E. NEWTON, SUMMER 1986 



V - MODAL MATRIX, N X N. (COLUMNS ARE EIGENVECTORS.: 
B - SUBDIAGONAL OF STIFFNESS MATRIX 
C - DIAGONAL OF STIFFNESS MATRIX 
D - DIAGONAL OF MASS MATRIX 

E - SUBDIAGONAL OF MASS MATRIX 

X - VECTOR OF EIGENVALUES 

N - VECTOR SIZE (=NI ) 



INTEGER J,K,N,L,N1 
PARAMETER(N1=12) 

REAL V(N,N),P(N1),Q(N1),T(N1),X(N),B(N),C(N),D(N),E(N), 
1 H, VN, TN, DIAG, XI , X2 , VTMV , TTMV , TTMT , SUM 



L=0 

DO 30 K=1,N 

IF(L. NE. 0)THEN 
L=0 

GO TO 30 
END IF 
X1=X(K) 

V( 1, K) =1 . DO 
V(2, K)=0. DO 
T(1)=0.D0 
T(2)=1.D0 
P(l)=B(l)-Xli^E(l) 

P(2)=B(2)-X1*E(2) 

VN=-(C( 1 )-Xl*D( 1 ) )/P( 1) 

TN=-P(2)/P(1) 

DO 10 J=2,N-1 

DI AG=C( J ) -X1»D( J ) 

P(J+1)=B(J+1)-X1»E(J+1) 

V( J+1,K)=-(P(J)»V( J-1,K)+DIAG»V(J,K) )/P( J+1) 
T(J+1)=-(P(J)»T(J-1)+DIAG»T(J))/P(J+1) 

10 CONTINUE 
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C CHECK FOR A DOUBLE ROOT 



1F(K.EQ.N)G0 TO 101 
X2=X(K+1) 

CRIT=ABS( (X2-X1)/(X2+X1) ) 

1F(CRIT.LT.4.E-4)G0 TO 22 

C CONSTRUCT EIGENVECTOR FOR SINGLE ROOT 

101 H=-(V(N,K)-VN)/(T(N)-TN) 

DO 11 J=2,N 

V(J,K)=V(J,K)+H»T(J) 

11 CONTINUE 

C NORMALIZE WITH RESPECT TO MASS MATRIX 

P(l)=D(l)«V(l,K)+E(2)»V(2,K)+E(l)i^V(N,K) 

DO 12 J=2,N-1 

P( J )=E( J ) J^V( J-1, K) +D( J ) »V( J , K) +E( J + 1 ) ifV( J + 1 , K) 

12 CONTINUE 

P(N) =E(N) i^V(N-l , K) +D(N) ifV(N, K) +E( 1 ) 1 , K) 

VTMV=0. DO 
DO 14 J=1,N 

VTMV = VTMV + P( J ) *V( J , K) 

14 CONTINUE 

VTMV=1 . DO/SQRT( VTMV) 

DO 16 J=1,N 

V(J,K)=V(J,K)»VTMV 
16 CONTINUE 

GO TO 30 

C CONSTRUCT EIGENVECTORS FOR DOUBLE ROOT AND NORMALIZE 

22 L=1 

P( 1 )=D( 1 )*V( 1 , K) +E(2) *V(2, K) +E( 1) »V(N, K) 

Q( 1 )=D( 1 ) »T( 1 ) +E(2) »T(2) +E( 1 ) »T(N) 

DO 23 J=2,N-1 

P(J)=E(J)»V(J-l,K)+D(J)»V(J,K)+E(J + l)ifV(J + l,K) 
Q( J ) =E( J ) i^T( J-1 ) +D( J ) i^T( J ) +E( J + 1 ) i^T(J + l ) 

23 CONTINUE 

P(N)=E(N)»V(N-l,K)+D(N)J(V(N,K)+E(l)<fV(l,K) 

Q(N)=E(N) »T(N-1 ) +D(N) *T(N) +E( 1 ) »T( 1 ) 

VTMV=0. DO 
TTMV=0. DO 
TTMT=0. DO 
DO 24 J=1,N 

VTMV=VTMV+P(J ) »V( J , K) 

TTMV = TTMV + P( J ) i^T( J ) 

TTMT=TTMT+Q( J ) »T( J ) 

24 CONTINUE 
H=-TTMV/VTMV 

SUM=TTMT+2. DO»H»TTMV+H»H» VTMV 
VTMV=1. DO/SQRT(VTMV) 

SUM=1 . DOZSQRT(SUM) 

H=H*SUM 
DO 26 J=1,N 

V(J,K + i:)=V(J,K)*H + T(J) i^SUM 
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V(J,K)=V(J,K)»VTMV 
26 CONTINUE 

30 CONTINUE 
RETURN 
END 



cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE RAYLYP ( X , V , B , C, P , N ) 



* SUBROUTINE RAYLYP USES THE RAYLEIGH QUOTIENT TO FIND IMPROVED 

* EIGENVALUES FROM THE ALREADY NORMALIZED EIGENVECTORS 






C 

C AUTHOR 
C 

C ARGUMENTS 



C OUT - X 

C IN - V 

C 

C B 

C C 

C P 

C N 

C 

C 



R. E. NEWTON, SUMMER 1986 



- VECTOR OF N EIGENVALUES IN NONDESCENDING ORDER 

- MODAL MATRIX, N X N. (NORMALIZED WITH RESPECT TO 
MASS MATRIX. ) 

- SUBDIAGONAL OF STIFFNESS MATRIX 

- DIAGONAL OF STIFFNESS MATRIX 

- WORK VECTOR 

- VECTOR SIZE (= NI ) 



INTEGER J,K,N 

REAL V(N, N) , P (N) , X(N) , B(N) , C(N) , XI 
DO 20 K=1,N 

P(1)=C(1)»V(1,K)+B(2)*V(2,K)+B(1)*V(N,K) 

DO 12 J=2,N-1 

P(J)=B(J)»V(J-1,K)+C(J)»V(J,K)+B(J+1)*V(J+1,K) 

12 CONTINUE 

P(N)=B(N)ifV(N-l,K)+C(N)»V(N,K)+B(l)*V(l,K) 

X1=0. DO 
DO 16 J=1,N 

X1=X1+P( J) *V( J,K) 

16 CONTINUE 

X(K)=X1 

20 CONTINUE 
RETURN 
END 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE SLVBVP(PHI , RHS , E I GVEC , B I GE, B I GC , B I GA , WK , D I R, HELM, 

1 N I , NJ ) 

» SUBROUTINE IS A SUBSTITUTE FOR EBVP2D 

C 

C AUTHOR - R. E. NEWTON, SUMMER 1986 

C 

C ARGUMENTS • 
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c 


OUT - PHI 


- 


SOLUTION 


c 


IN - RHS 


- 


RIGHT-HAND SIDE 


c 


El GVEC 


- 


MATRIX OF EIGENVECTORS 


c 


B IGE 


- 


INVERSE OF DIAGONAL FACTOR OF COEFFICIENT MATRIX 


c 

c 


B IGC 




SUPERDIAGONAL OF UNIT UPPER TRIANGULAR FACTOR OF 
COEFFICIENT MATRIX 


c 

c 


B IGA 




SUBDIAGONAL OF UNIT LOWER TRIANGULAR FACTOR OF 
COEFFICIENT MATRIX 


c 


WK 


- 


WORK AREA 


c 


DIR 


- 


LOGICAL SWITCH - TRUE FOR DIRICHLET B.C. 


c 


HELM 


- 


LOGICAL SWITCH - TRUE FOR HELMHOLTZ EQUATION 


c 


NI 


- 


X-DIMENS ION 


c 

r 


NJ 


- 


Y-DIMENSION 



REAL PHI (NI , NJ ) , RHS(NI , NJ ) , EIGVEC(NI , N I ) , B I GE ( N I , N J ) , B 1 GC ( N I , NJ ) , 
1 B IGA(NI , NJ ) , WK(NJ ,NI) , DUM 
LOGICAL DIR, HELM 

C TRANSFORM RIGHT-HAND SIDE 



CALL MATPR ( WK, RHS, E I GVEC, NJ , N I , N I , NJ , N I , N I , DUM, DUM, DUM, D I R, 
1 . FALSE. ) 

C PERFORM FORWARD REDUCTIONS 



5 



10 

15 



DO 5 I = 1, NI 

WK(1,I) = WK(1, I )*BIGE( I , 1) 

CONTINUE 
.DO 15 I = 1, NI 

DO 10 J = 2, NJ 

WK(J,I) = WK( J , I ) »B IGE( I , J ) + WK(J-1 
CONTINUE 
CONTINUE 



»B IGA( I , J ) 



C PERFORM BACK SUBSTITUTIONS 

DO 25 I = 1, N I 

DO 20 J = NJ-1, 1, -1 

WK(J,I) = WK(J,I) + WK(J+1, I )»BIGC( I , J) 
20 CONTINUE 

25 CONTINUE 



C BACK TRANSFORM RESULTS 

DO 40 J = 1, NJ 

DO 35 I = 1, NI 
DUM = 0. 

DO 30 K = 1, NI 

DUM = DUM + EIGVECC I ,K)»WK( J,K) 
30 CONTINUE 

PHI ( I , J ) = DUM 
35 CONTINUE 

40 CONTINUE 
RETURN 
END 
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cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE TRIEIG(V,X,B,C,D,E,WU,N) 

» SUBROUTINE TR I E I G USES THE METHOD OF BISECTION TO FIND EIGENVALUES » 

# FOR THE GENERALIZED EIGENVALUE PROBLEM INVOLVING SYMMETRIC TRIDIA- 

# GONAL MATRICES. THE ROUTINE IS ADAPTED FROM AN ALGOL PROGRAM NAMED# 
» 'BISECT' WRITTEN BY BARTH, MARTIN, AND WILKINSON (NUM. MATH. 9, 386- 

# 393 (1967)). 'BISECT' FINDS EIGENVALUES (STANDARD EIGENVALUE PROB- 

# LEM) FOR A SYMMETRIC TRIDIAGONAL MATRIX. THE MATRICES FOR TR I E I G 

# ARE INPUT AS VECTORS - C(N) AS THE DIAGONAL AND B(N) AS THE SUB- « 

# DIAGONAL OF THE "STIFFNESS" MATRIX AND D(N) AS THE DIAGONAL AND E(N) 

# AS THE SUBDIAGONAL ON THE "MASS" MATRIX. SUBROUTINE EIGVEC FINDS 

# THE EIGENVECTORS AND NORMALIZES THEM WITH RESPECT TO THE MASS MATRIX 

c 

C AUTHOR: 

C 

C ARGUMENTS 
C OUT - 

C 
C 

C IN - 

C 
C 
C 
C 
C 
C 

c 

c 

INTEGER N, Ml , N, Z, I , A, K 

REAL C(N),B(N),X(N),WU(N), EPSl , EPS2, RF, FI , XMI N, XMAX, 

IQ, XI , XU, XO, DELT, D(N) , E(N) , DC, DB , DD , DE, TMAX, TMIN,V(N,N),G(N,N) 
DATA EPSl , RF/0. , 1 . E-7/ 

C CALCULATION OF XMAX, XM I N 

B( 1 ) =0. 

E( 1 ) =0. 

DC=C(N) 

DB=ABS(B(N) ) 

DD=D(N) 

DE=E(N) 

XMAX=(DC+DB)/(DD+DE) 

XMIN= (DC-DB) / (DD-DE) 

C NEGATIVE EIGENVALUES ASSUMED. IF EIGENVALUES ARE ALL POSITIVE, 

C REPLACE THE TWO PRECEDING LINES WITH THE FOLLOWING TWO LINES. 

# XMAX= (DC+DB) / (DD-DE) 

» XMIN=(DC-DB) / (DD-t-DE) 

DO 2 I=N-1, 1, -1 
DC=C( I ) 

DB=ABS(B( I ) )+ABS(B( I+l) ) 

DD=D( I ) 

DE=E( I ) +E( I +1 ) 

TMAX=(DC+DB) /(DD+DE) 



R. E. NEWTON, SUMMER 1985 



V 

X 

C 

B 

D 

E 

WU 

N 



MATRIX OF EIGENVECTORS, N X N, NORMALIZED WITH 
RESPECT TO THE "MASS" MATRIX 

VECTOR OF EIGENVALUES IN NONDESCENDING ORDER 
DIAGONAL OF "STIFFNESS" MATRIX 

OF "STIFFNESS" MATRIX 
"MASS" MATRIX 
OF "MASS" MATRIX 



SUBDIAGONAL 
DIAGONAL OF 
SUBDIAGONAL 
WORK VECTOR 
VECTOR SIZE 



( = N I ) 
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TMIN= (DC-DB) / (DD-DE) 

C NEGATIVE 'eigenvalues ASSUMED. IF EIGENVALUES ARE ALL POSITIVE, 
C REPLACE THE TWO PRECEDING LINES WITH THE FOLLOWING TWO LINES. 

* TMAX=(DC+DB)/ (DD-DE) 

» TMIN= (DC-DB) / (DD+DE) 

IF(TMAX. GT.XMAX)XMAX=TMAX 
I F ( TM I N . LT . XM I N ) XM I N=TM I N 
2 CONTINUE 

C SET EPS2 

IF(XMIN+XMAX.GT.O. ) THEN 
EPS2=RF»XMAX 

ELSE 

EPS2=RF* (-XMIN) 

END IF 

IF(EPSl.LE.O. )EPS1=EPS2 
EPS2=0. 5»EPSl+7. »EPS2 

C INNER BLOCK 

XO=XMAX 
DO 4 I =N, 1 , -1 
X( I )=XMAX 
WU( I )=XMIN 
4 CONTINUE 

C LOOP FOR K-TH EIGENVALUE 

DO 100 K=N, 1,-1 • 

XU=XMIN 

DO 6 I =K, 1 , -1 

IF(XU.LT. WU( I ) )THEN 
XU=WU( I ) 

GO TO 10 
END IF 
CONTINUE 

IF(XO.GT.X(K) )XO=X(K) 

Xl= (XU+XO) /2. 

Z = Z+1 

C STURM’S SEQUENCE 



6 

10 

20 



A = 0 
Q=l. 

DO 30 1=1, N 

IF(Q.EQ. 0. )THEN 

DELT=ABS(B( I )-Xl»E( I ) ) /RF 

ELSE 

DELT=(B( I )-Xl*E( I ) )»*2/Q 
END IF 

Q=C( I ) -X1»D( I ) -DELT 
IF(Q. LT.O. )A=A+1 
30 CONTINUE 

IF(A. LT. K)THEN 

IF(A. LT. 1 )THEN 
WU(1)=X1 
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XU = X1 



ELSE 

WU(A+1 ) =X1 
XU=X1 

IF(X(A) . GT. XI )X(A) =X1 
END IF 

ELSE 

X0 = X1 

END IF 

IF(X0-XU.GT.2. *RF^KABS(XU)+ABS(XO) )+EPSl)GO TO 20 
X(K)=(X0+XU)/2. 

100 CONTINUE 

CALL EIGVEC(V,B,C,D,E,X,N) 

CALL RAYLEE(V,B,C,X,N) 

CALL E IGVEC(V, B, C, D, E, X, N) 

RETURN 

END 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX> 
SUBROUTINE E I GVEC ( V , B , C , D , E, X, N ) 

* SUBROUTINE EIGVEC FINDS EIGENVECTORS BY DIRECT SOLUTION OF THE GOV- 

* ERNING LINEAR ALGEBRAIC EQUATIONS AND NORMALIZES THEM WITH RESPECT 
» TO THE MASS MATRIX. 

C 

C AUTHOR 
C 

C ARGUMENTS 
C OUT - 

C IN - 

C 
C 
C 
C 
C 

C 

C 

INTEGER J,K,N 

REAL V(N,N),P(5),X(N),B(N),C(N),D(N),E(N),X1, DSQRT, SUM 

DO 20 K=1,N 
X1=X(K) 

V(1,K)=1. 

P(2)=B(2) -E(2) »X1 
V(2,K)=(D(1)*X1-C(1))/P(2) 

DO 10 J=2,N-1 

P(J+1)=B(J+1)-E(J+1)*X1 

V(J+1,K)=-(P(J)#V(J-1,K)+(C(J)-D(J)*X1)*V(J,K))/P(J+1) 

10 CONTINUE 

C NORMALIZE WITH RESPECT TO MASS MATRIX 

P(1)=D(1)*V(1,K)+E(2)*V(2,K) 

DO 12 J=2,N-1 



- R. E. NEWTON, SUMMER 1985 



V - MODAL MATRIX, N X N. (COLUMNS ARE EIGENVECTORS. 
B - SUBDIAGONAL OF STIFFNESS MATRIX 
C - DIAGONAL OF STIFFNESS MATRIX 

D - DIAGONAL OF MASS MATRIX 

E - SUBDIAGONAL OF MASS MATRIX 

N - VECTOR SIZE (= NI) 
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P(J)=E(J)»V(J-1,K)+D(J)»V(J,K)+E(J + 1)»^V(J + 1,K) 
CONTINUE 

P(N)=E(N)»V(N-1,K)+D(N)*V(N,K) 

SUM=0. 

DO 14 J=1,N 

SUM=SUM+P(J)*V(J,K) 

14 CONTINUE 

SUM=1. /SQRT(SUM) 

DO 16 J=1,N 

V(J,K)=V(J,K)»SUM 
16 CONTINUE 

20 CONTINUE 
RETURN 
END 



cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE RAYLEE-( V, B, C, X, N) 



****1i***H*************it***************1Hi****************^*^************* 

* SUBROUTINE RAYLEE USES THE RAYLIEGH QUOTIENT TO FIND IMPROVED EIGEN-* 

* VALUES FROM THE ALREADY NORMALIZED EIGENVECTORS. * 

it*********************************************************************** 



C 

C AUTHOR - R 

C 

C ARGUMENTS 
C OUT - X 

C IN - V 

C 

C B 

C C 

C N 

C 

C 



E. NEWTON, SUMMER 1985 



VECTOR OF EIGENVALUES IN ASCENDING ORDER 

MODAL MATRIX, N X N. (NORMALIZED WITH RESPECT TO 

MASS MATRIX) 

SUBDIAGONAL OF STIFFNESS MATRIX 
DIAGONAL OF STIFFNESS MATRIX 
VECTOR SIZE ( = NI) 



12 



16 

20 



INTEGER 


J , K, N 










REAL V(N 


,N),P(5), 


X(N) , 


B(N) 


, C(N) 


,X1 


DO 20 K= 


1,N 










P(1 ) = 


C(1)*V(1, 


K) +B( 


2) *V 


(2,K) 




DO 12 


J=2,N-1 












P(J)=B(J) 


*V( J- 


1,K) 


+ C( J ) 


*V(J,K)+B(J+1)*V(J+1,K) 


CONTI 


NUE 










P(N) = 


B (N) *V(N- 


1,K) + 


C(N) 


*V(N, 


K) 


X1=0. 












DO 16 


J = 1,N 












X1=X1+P( J 


) »V( J 


,K) 







CONTINUE 

X(K)=X1 

CONTINUE 

RETURN 

END 
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